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Abstract In the weak field limit general relativity reduces, as is well known, to the 
Newtonian gravitation. Alternative theories of gravity, however, do not necessarily 
reduce to Newtonian gravitation; some of them, for example, reduce to Yukawa-like 
potentials instead of the Newtonian potential. Since the Newtonian gravitation is 
largely used to model with success the structures of the universe, such as for example 
galaxies and clusters of galaxies, a way to probe and constrain alternative theories, in 
the weak field limit, is to apply them to model the structures of the universe. In the 
present study we consider how to probe Yukawa-like potentials using A^-body numerical 
simulations. 

PACS 04.50.+h, 04.50.Kd, 45.50.-j 



1 Introduction 



Einstein's General Relativity (GR) is one of the most beautiful theories ever imagined 
by the human mind. Although GR is a successful theory of gravitation, it is unable to 
explain, for example, the accelerating expansion of the universe, unless a cosmological 
constant or a dark energy fluid is considered. Nowadays other theories of gravitation 
intend to give an alternative interpretation to this accelerating expansion (see, e.g., 
PQ). Most of these alternative theories, although having different approaches (some 
scalar-tensor theories of gravitation, nonsymmetric gravitational theory, etc), reduce, 
in the weak-field limit, to a Yukawa-like gravitational potential (hereafter YGP), i.e., 

0( r ) = _GM_ e - r /x^ (1) 
r 
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The above equation gives the potential of a point mass m at a distance r; G is the 
universal gravitational constant and A is a constant. When A — » oo, this potential tends 
to the Newtonian one. Note that the parameter A is the Compton wavelength of the 
exchange particle, which in the present case is a graviton. The graviton mass is related 
to A through the well known equation m g = h/Xc, where h is the Planck constant and 
c is the speed of light. 

The YGP has been investigated in the literature, in particular, in galactic 
astronomy and cosmology. For example, Signore [2] studies this potential under 
the cosmological context, on the other hand, de Araujo & Miranda [3] study how 
variations of the A parameter can disturb galactic disks. This last study generalizes 
the (Newtonian) Toomre's [4] expressions for rotation curves in order to account for a 
YGP. Rodriguez-Meza et al. [5] consider the YGP as a correction on the Newtonian 
potential. 

Recently, some numerical approaches have been developed to study alternative 
theories. Cervantes-Cota et al. [6j[7] have developed numerical studies of a scalar-tensor 
theory in the weak field limit, using the iV-body method. In these studies they consider 
that the gravitational potential in the weak field limit is given by 

&STT = <t>N(f) + acj){r)/{\ + a) , 

where $stt is the potential from the scalar-tensor theory, <f>N( r ) is the Newtonian 
potential, 4>(r-) is given by Eq.(l), and a = 1/(3 + 2ui), with uj being the Brans- 
Dicke parameter [8] . They changed a particular iV-body code replacing the Newtonian 
potential by the potential <&stt an d analyzed its influence on isolated galaxies, two 
colliding spiral galaxies and issues concerning the formation of bars. Simulations of 
cosmological structure formation in ACDM scenarios are also studied. 

In the present paper, we use similar techniques to that used by Cervantes et al. [6| 
[7], although with a gravitational potential derived from other theories. We modify a 
popular iV-body code, Gadget-2 9 , replacing the Newtonian potential by a pure YGP 
as given by Equation [T] which is derived from some alternative theories of gravitation 
(see, e.g., Visser's theory [lOjV 

It is worth stressing that although we have chosen a particular gravitational 
potential, the approach considered here can be applied to any alternative gravitational 
potential. 

This paper is organized as follows: in section 2, we show the changes made in the 
Gadget-2 code, in section 3, we test the YGP Gadget-2 code using a pair of particles 
and iV-body systems. Finally, in section 4, we present the conclusions and briefly discuss 
how this modified code can be used to test alternative theories of gravitation via galactic 
dynamics. 

2 Changing the Gadget-2 Code 

Gadget-2 code is a new massively parallel TreeSPH Code, designed to simulate the 
evolution of iV-body collisionless systems, like, e.g., cosmological fluids, interacting 
galaxies, etc. It was developed to work with the concept of particles: matter distribution 
is represented by a system of particles, where each particle has position, velocity, label, 
gravitational potential, density (if it represents a gas), and so on. 

Gadget-2 uses an advanced system of processor's communication, based on MPI 
interface. MPI is a computer implementation that allows many computers to work in 



3 



parallel. Thus, one can simulate from thousand of particles using desktops machines 
to billions of particles using supercomputers. This code computes quickly many 
gravitational forces using a hierarchical tree algorithm, like the Barnes' Tree Code 
[11] . grouping distant particles into larger cells and computing their total potential 
by multipole expansions. This method is very efficient and the computational effort 
is 0(N log N) interactions, faster than the N 2 interactions from direct sum method. 
Gadget-2 also computes gravitational potentials using particle-mesh (PM) method, 
when cosmological simulations are considered. In this case, the so called "long range" 
potentials are calculated by the PM method, which computes potentials in the Fourier 
space using the mesh method. On the other hand, the "short range" potentials are 
calculated in real space using the tree algorithm. 

It is worth stressing another very important characteristic of the Gadget-2 code: 
its capability to integrate the equations of motion of particles with a leapfrog integrator 
based on the the KDK scheme (KDK means kick-drift-kick, see [9]), which conserves the 
symplectic nature of collisionless systems, in other words, it preserves the Hamiltonian 
structure of these simulated systems. This approach is more efficient than other ones 
if applied with the adaptive timestep method, as shown by Springel [9]. In this way, 
Gadget-2 is in general more efficient than codes based on the Runge-Kutta integrators. 

It is worth mentioning that Gadget-2 offers a set of options: pure newtonian 
simulations, cosmological comoving periodic simulation with the alternative use of 
the TreePM formalism, etc. However, the "Cosmological simulation options" can be 
activated or not, whether the code is compiled with these respective options. To 
simulate only gravitational interactions in local systems, such as galaxies or globular 
clusters, we compile this code without the advanced cosmological options listed in the 
Makefile, "turning on" the gravitation, and "turning off' all the advanced options, such 
as PM formalism, expansion of the universe, etc. 

In the present paper, we intend to study only gravitational interactions in local 
systems, without gas and cosmological expansion. To this end, we compile Gadget-2 
without the advanced options, and so we are able to study a cloud of particles self- 
interacting gravitationally. 

We fully studied the Gadget-2 code in order to see where it is necessary to modify 
it in order to take into account the corresponding YGP expressions. We concluded that 
it is only necessary to change the code where instructions related to the Newtonian 
potential and accelerations appear. In this case, the acceleration reads 

a = -^e-A(r + A )r . (2) 
It is important to note that we have left A as a free parameter. 



3 Testing the Code 

The reliability of a JV-body simulation to model a given system depends on how one 
can control the numerical errors naturally intrinsic in this approach. 

For example, the particles have their accelerations computed by a tree recursive 
code and, consequently, their positions and velocities are updated from the integrators. 
As a result the total energy of the simulated system, can be well conserved or not, 
whether the treecode calculations and the integration processes are well accurate or 
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not. For these reasons, errors must be investigated and minimized to an acceptable 
level. 

In the TV-body numerical formalism, one uses the fact that gravitation obeys the 
superposition principle. Therefore, one should firstly test the two-body problem, testing 
the code's integrator precision to the new gravitational potential one is considering. 
Subsequent tests should consider TV-body systems, verifying not only the treecode 
precision, but also the influence of some numerical parameters on the timestep values. 

3.1 Simulating the motion of two particles 

In this first test, we are able to probe how our changes in the Gadget-2 code modify 
its integrator's precision. A similar test was performed by [9] using newtonian physics. 
This test is very important, because we are interested in verifying the code's precision 
to drift positions and velocities in each timestep for the YGP. 

On the other hand, one could argue that TV-body systems and a two-body system 
are very different numerically, and consequently that this test would be useless. 
However, the TV-body formalism is based on the superposition principle, in other words, 
the total potential is the sum of the potential from all individual particles. Then, we 
firstly test the two-body problem in order to verify if the code can work appropriately. 
Subsequently, we performe tests with a large amount of particles. 

In fact, in the case where treecodes are used, the superposition principle is used 
implicitly in the monopole, quadrupole and octopole expansion terms. But, Gadget- 
2 uses only the monopolar expansions, which is sufficient and very accurate in the 
newtonian case, as shown by [9]. In the case where the YGP is used, we have performed 
some tests and, as we will show in the next sections, its monopolar term is also enough 
to yield satisfactory results. 

In this first test, we simulate the motion of two particles: a massive particle and a 
test particle (whose mass can be neglected). The massive particle has mass M=100.0 
(in code units, i.e., M=1.0 = 1.0 x 1O 1O M0), and is located at (0,0,0) kpc, the origin 
of a cartesian coordinate system, and has velocity v = (0, 0, 0) km.s -1 . It is important 
to note that this system of units is very common in galactic simulations with the 
Gadget-2 code. In this system, the time unit is given in 9.8 x 10 s ~ 10 9 years. 

The test particle has mass m = 0.1, and is located initially at r = (100,0,0) kpc 
with v p = (0,104,0) km.s" 1 . We have chosen A = 100 kpc, for our computational 
purposes. Note that the initial separation between the particles is equal to A. In this 
simulation, the timestep At values are 0.00146484 < At < 0.00585938. It is worth 
noting that in Gadget-2, the timestep criterion is 



where n is an accuracy parameter; Atmax, the the maximum timestep allowed, which 
is usually chosen to be a small fraction of the dynamical time of the system under 
study; e, the gravitational softening; and ja|, the modulus of the particle's acceleration, 
estimated from the last timestep. In the subsection 3.1.1, we consider the role of the 
timestep in our simulations. 

We have made the system orbiting for 6 Gyears. So, we are able to study the 
trajectory of the test particle in detail. We display the test particle's orbit in Fig. [T] 
This picture is made by displaying the particle's position in the xy— plane obtained 
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Fig. 1 Trajectories of the particles on the xy— plane. The central particle is the massive one. 



from 600 snapshots, including the initial one. As can be seen, the test particle orbits 
in the xy— plane and precesses. This orbit characteristic is very important to test the 
reliability of our modified code. 

The orbit precession showed in Fig.[l]is a real physical effect, it is not related to any 
numerical instabilities. YGP decreases with distance more steeply than the Newtonian 
potential. It is well known, from the Bertrand's theorem, that potentials that fall with 
distance differently from 1/r make the orbits to precess; in other words, the orbits are 
not closed, as we show in Fig. [T] In the following sections, we will study this issue in 
detail. 



3.1.1 Timesteps in the two-body simulation 

From Equation [3] the At values depend strongly on the acceleration of the particles. 
But the acceleration depends also on some other parameters, for example, the softening 
length e and the YGP parameter A. When the distance between two particles is smaller 
than the softening length scale, the gravity law in the simulation is replaced by a 
softening kernel that avoids divergences in the calculated field (see, e.g., [5]). Due to 
the numerical nature of our calculations, the investigation of the timestep values is a 
very important task, due to the fact that their values are used by the KDK algorithm 
to update positions and velocities. Therefore, it is necessary to probe the influence 
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Table 1 Some runs with different e and A values. The first column shows the A values, the 
second one the e values and the third one, At. The first line display the results of the simulation 
described in Subsection 13. II 
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1.0 


0.0059 


1.0 


0.1 


0.0059 



of the numerical parameters A and e on the values of the timesteps. In this sense, 
we have rerun the two-body simulations described above, changing the values of the 
parameters A and e. Concerning the iV-body system, we will show in the next sections 
why Gadget-2 is a safe code to deal with these parameters, displaying the At values 
for different parameters A and e. 

In Table [1] we display the maximum At values from a set of 6 runs, where we 
have set the same simulated time and used the same two-body system described in 
Subsection 13.11 We have only changed the values of e and A, setting n — 0.025 for all 
experiments, to improve the timestep accuracy. 

In these simulations, we run the two-body system for 6 Gyears for A = 1, 10 and 
100 kpc, and with e — 0.1 and 1 kpc. Although setting e values for two-body tests 
seems to be unnecessary, we want to probe how these parameters change our results. 
We have also made a newtonian simulation for a two-body system using the default 
Gadget-2, to compare our modified code to the original one. In the newtonian case, 
we obtain 0.0015 < At < 0.0031 for e = 1 kpc and 0.00038 < At < 0.00076 for e = 0.1 
kpc. 

It is worth mentioning that for all the runs, for different values of A, e and timesteps, 
the total energy conserves very well, as we will see below. 

3.1.2 Testing the reliability of our simulations 

A good integrator must conserve the total energy of the simulated system within 
an acceptable level. However, due to the nature of the YGP, orbits are not closed 
and this feature reveals itself as a good indicator to probe how reliable is our two- 
body simulation. So, we performed four additional tests to verify the reliability of 
our simulation described in the Subsection 13.11 namely: energy conservation, angular 
momentum, orbital period and precession of the orbit. 

In the energy conservation test, we compare our results with the Newtonian case. 
Thus, we have made another simulation with the original (Newtonian) code. This 
simulation is performed with the same initial conditions, as we explained in the last 
section. In Figs. [2] and [3] we show the total energy evolution for both simulations. 
These pictures show the evolution of the total energy E normalized by the total initial 
energy Eg. To the cases where A = 100 kpc, we consider e = 0.1 and e = 1 and the 
results are practically the same. Therefore, we display here only the results for e = 1 
and A = 100 kpc. 

In Fig. [2] we conclude that the maximum energy violation (E — Eq)/Eq is about 
0.6% at t ~ 3.3 Gyears. In the YGP case, as can be seen in Fig. [3] the energy violation 
is systematically (E — Eq)/Eq ~ 1.3%. Note that the energy conservation is better in 
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Fig. 2 The total energy (normalized by the total initial energy Eq) evolution as function of 
time for the Newtonian case. 



Newtonian case, this has to due with the fact that 1/r behaves numerically better than 
g -r/A jj oweverj we consider that this violation can be neglected for two reasons: firstly, 
6 Gyears is a great time lapse for such an energy violation. Secondly, this violation is 
systematic, at some snapshots the violation is null and this behavior also appears in 
Newtonian case. It is interesting to note that at 6 Gyears the energy violation is almost 
null in both situations, showing that our method is consistent. 

One could ask what happen with the energy conservation if we follow the 
simulations, say, till 100 Gyears. The results of a such simulation show that the energy 
conservation is systematically small, just like in the 6 Gyears simulation. Such a result 
shows again the reliability of our code. 

We have also tested the conservation of angular momentum. Fig. U shows the 
evolution of the total angular momentum L normalized by the total initial angular 
momentum Lq. We note that the total angular momentum is violated by an amount 
of 0.2 % over the initial value. 

Besides the energy and the angular momentum conservation, we have also 
performed three additional tests. The first one is to verify the orbital half-period of 
the test particle. The second one is to check the precession of the particle's orbit. The 
third test investigates the difference between the potentials calculated by the tree and 
the direct sum methods in the modified code. To do this test, we use 30,000 particles 
to represent a typical snapshot of galaxy simulations. 
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Fig. 3 The total energy (normalized by the total initial energy Eq) evolution as function of 
time for the Yukawian case with A = 100 kpc. 



In the test of the period, we evaluate the half-period, which is given by: 



where T r is the orbital half-period, r a p is the apocenter, r p is the pericenter, E is the 
initial total energy per unit mass, <2>(r) is the potential over the test particle, r is the 
distance from the center of mass to the test particle and L the angular momentum 
per unit mass. We calculate the above integral numerically and found a value of 0.78 
Gyears. In our simulation, the numerical value obtained from the snapshots is also 0.78 
Gyears, i.e., the results agree quite well. This result can also be seen in Figs. [5] and [6] 
The second test consists in the evaluation of the precession angle. Bertrand's 
theorem states that potentials decreasing with distance differently from 1/r make the 
orbits to precess (see, e.g., Goldstein |12|). We then evaluated the precession using the 
following equation: 



where Aip is the shift of the semi major axis. We again evaluate this integral numerically 
and found that Aip — 28.72°. This value agrees with the results obtained from a 
snapshot of our simulation. Taking the snapshot at t=0.78 Gyears, we found that the 
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Fig. 4 Evolution of the angular momentum (normalized by the initial angular momentum) 
for A = 100 kpc simulation. 

particle's position, from the first snapshot, is A$ — 28.72°. This result again shows 
that our YGP Gadget-2 works quite well and reproduces many physical properties of 
this kind of potential, due to the good accuracy provided by its integrator's algorithm. 

3.2 Tests with iV-body systems 

So far we have investigated the two-body problem, whose tests demonstrated that our 
modified code produces reliable results. Now, we consider the iV-body problem, for 
which Gadget-2 was designed. In this way, we are initially interested in the recursive 
octaltree algorithm accuracy to compute the YGP, when compared to the particle- 
particle method. We are also interested in the influence of the numerical parameters 
N, A and e over the timestep values, as well as in the energy conservation. In what 
follows, we will show the main results of a set of iV-body simulations. 

3.2.1 The recursive tree algorithm in the YGP case 

These tests consist in the comparison of the gravitational potentials of a system of 
particles calculated by two different methods: the exact one, also called the direct sum 
method or particle-particle (PP) method, that computes all potentials using ~ TV 2 
operations, and the treecode method, used by Gadget-2 code, as we explained above. 
In this test, we use both codes, the default Gadget-2 and our modified version. To 
do this, we have chosen a typical snapshot used in galactic dynamics. This snapshot 
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Fig. 5 Evolution of the particle's distance from the center of mass frame for the A = 100 kpc 
simulation. 



describes a system of particles modeling the initial conditions of one exponential 
disk galaxy composed by 10,000 particles embedded in a Hernquist dark matter halo 
composed by 20,000 particles. Both systems are built from typical parameters, like the 
galaxy models developed by Springel & White |13| . Later on we will see that such a 
resolution is enough to test the reliability of the modified treecode method to compute 
the YGPs of all the particles. 

In Fig. [7] we show a visualization of our galaxy model. In another publication to 
appear elsewhere we show such galaxy models in detail. 

Using a Fortran code, we have computed the total potential of each halo and 
disk particles by the direct sum method and compared these results to the potentials 
computed by Gadget-2 code. We define A<f> = \4>Nbody ~ 4*Tree\ as the modulus 
of the difference resulting from computation of the potential by the two different 
methods. 4>Nbody indicates the potential calculated by the direct sum method with our 
Fortran code, and 4>Tree, the potential calculated with the Gadget-2, that uses the 
tree method. To do this, we have used typical parameters of the tree, e.g., the tolerance 
parameter 6 — 0.8, and softening lengths of e/, = 1 kpc for the dark matter halo and 
e<j = 0.4 for the baryonic disk. These parameters are unchanged in the Newtonian and 
the YGP cases, i.e., they are fixed parameters. Defining the relative error A<f>/<f>Nbodyi 
we can compare both methods. The system of units used here is 1O 1O M0 for unit mass, 
kpc for distances and km.s^ 1 for velocities. 
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Fig. 7 Typical snapshot modeling a spiral galaxy and its dark matter halo seen in perspective. 
Each box side has 80 kpc. Blue particles represent dark matter halo and white/yellow particles 
indicate the disk. 
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Fig. 8 Modulus of the relative errors as a function of the gravitational newtonian potential 
for the particles of the dark matter halo and for the baryonic disk. 

In Fig. [8] we display A4>/4>Nbody against the potentials of the particles in the 
newtonian case, while in Fig.[9]we show the YGP case, having in mind that both cases 
are calculated to the same snapshot. To compute the potentials in the YGP case, we 
have set A=100 kpc, the same halo's length displayed in Fig. [7] Due to the fact that the 
halo and A have equal length in our snapshot, the e~ r / A factor plays a important role 
in the potential computations at medium and great separations between the particles. 
Note that if two particles are at distance of 2.5 kpc, the potential is 0.976 Newtonian 
(see Eq. [TJ. But, in the case where the distance is 10 kpc, the Yukawian potential is 
0.9 of the Newtonian, and when the distance is 25 kpc, the potential is 0.78 of the 
Newtonian. 

Analyzing Figs. [8] and [9] we conclude that in the YGP case for A = 100 kpc the tree 
precision behaves very similarly to the Newtonian case. Errors are the same in both 
cases and our figures show that the modified code can compute potentials of JV-body 
systems with the same precision as the Default Code. Note that we have used the initial 
snapshots in Figs. [8] and [9] but similar results are obtained for any other snapshots of 
our simulations. 

3.2.2 Additional N-body tests: A = 1 kpc and A = lOkpc. 

In addition to the tests described above, we make simulations now considering A = 1 
and 10 kpc, using a snapshot similar to the previous one, but with a better resolution, 
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Fig. 9 The same as Figure |8] for YGP with A = lOOkpc. 



namely: 30,000 particles for the disk and for the dark halo, totalizing 60,000 particles. 
In this new test, we have changed the e values to e/j = 0.1 and = 0.05. Our 
computational procedures are very similar to that described in the previous Subsection 
and here we show that 4>tree and (pNbody are practically the same. 

The results are displayed in the Figs. [TUJ and 1111 where we have plotted <f>tree 
against 4>Nbody f° r al1 the particles of our simulations; the dashed line represents 
4>Nbody = 4>tree- Note that the points are distributed along the dashed line, with a 
very small scattering. This result shows that our modified code calculates very well the 
gravitational potential. 



3.2.3 Timestep tests 

Numerical simulations are subject yet to considerable controversy, mostly concerning 
their numerical collisionless nature and their resolution. It is well known that the 
parameters N, e and At can change the simulation's results, and, in the galactic 
dynamic scenarios, the final morphology can depend on this set of parameters. 

An important detailed study in this issue is provided by Gabbasov et. al.[15|. where 
they probe the influence of the numerical parameters N, e and At on the outcomes of a 
set of simulations with disk galaxies and issues related to the formation of bars in these 
late-type systems. They prove, for example, that different sets of these parameters can 
change significantly the results of simulations that have initially the same snapshot. 
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Fig. 10 The octtree accuracy for A = 1 kpc case. (f>tree vs. (pNbody for all the particles of our 
simulations; the dashed line represents 4>Nbody = <frtree- 



In this paper, we make some tests similar to those by Gabbassov's team to probe 
the influence of the numerical parameters N, the total number of particles, e, and A on 
the timesteps values, using a methodology similar to that used in [15]. It is important 
to bear in mind that our tests are somewhat less complex and we performed a set of 
numerical simulations up to 0.1 Gyears of simulated time only to verify the influence 
of these parameters. 

This simulated time is short, if compared to the many runs performed by [15] . 
However, due to the Gadget-2 method in computing the timesteps, each of our runs 
take ~ 500-5000 steps to be completed, being a very satisfactory number of steps to 
our purposes. 

Besides, in some tests performed with 30,000 and 300,000 particles over 1 Gyear, 
the AE/Eq values did not change significantly during all the simulated time, although 
50,000 steps were necessary to this end. 

In this sense, 0.1 Gyears are sufficient for our studies. In this way, we have used three 
snapshots with 30,000, 300,000 and 1,200,000 particles, conserving the ratio 1:2 for the 
baryonic disk and dark halo, respectively. We run these snapshots for different values 
of e and A. These snapshots represent an exponential disk embedded in a Hernquist 
Sphere, as we explained above. 



15 




Fig. 11 The octtree accuracy for A = 10 kpc case, (fitree vs. <j>Nbody for all the particles of 
our simulations; the dashed line represents where (pNbody = </>tree- 

Table 2 Some runs for different N, e and A values (in kpc). Note that the respective At values 
are calculated by the Gadget-2 code. The energy conservation AE/Eo is also shown for each 
run. 
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0.004 


4.88281 xl0~ 5 


0.04% 


30 000 


100.0 


1.0 


0.4 


0.00078 


0.08% 


30 000 


100.0 


0.1 


0.04 


0.00020 


0.08% 


300 000 


100.0 


0.1 


0.04 


0.00020 


0.06% 


300 000 


100.0 


0.01 


0.004 


2.44141 xl0~ 5 


0.2% 


1 200 000 


100.0 


0.1 


0.04 


0.00020 


0.04% 


1 200 000 


100.0 


0.01 


0.004 


4.88281 xl0~ 5 


0.05% 
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To make the runs displayed in the Table El we estimate the e values taking the 
mean inter-particle separation at the center of the halo and disk. Our values are near 
to those used by [15] and are typical, if compared to those cited in the literature. Note 
that in the cases where e < 0.01, the simulations gets "more collisional" and the energy 
conservation gets a little bit worse, in agreement to the results obtained by [15]. If the 
resolution is increased, in other words, if the e values are smaller, then timestep values 
are smaller as well. As a consequence, the Gadget-2 tree algorithm computes greater 
accelerations and potentials at these smaller scales than in the small resolution cases 
(greater e values). 

However, in these cases where N is very large, the system tends to be 
collisionless, because in the close encounters, the softening kernel prevent divergences 
of accelerations and potentials, once e is chosen appropriately. Also, the collective 
potential dominates over individual interactions, making results better to N > 10 5 
particles. 

So, the greater the resolution is (i.e., the greater N) the better is the energy 
conservation as long as the e values are correctly chosen. But, note that, even with 
a modest value of N, the energy conservation can be quite good. Also, to probe a given 
alternative theory it is not always necessary a high resolution simulation, since the 
global properties of a galaxy, such as, for example, the density profile, can be very 
sensitive to the alternative potential. In other papers to appear elsewhere we consider 
this issue for early and late type galaxies. 

From our results, we consider our modified code able to test the YGP. Therefore, 
we are able to study alternative theories using galaxies or globular clusters, with the 
same quality one obtains with the use of Newtonian codes. 

4 Conclusions 

The present paper considers the use of the well known Gadget-2, with the appropriate 
modifications, to study alternative theories of gravitation, with particular attention to 
the YGP. In this first paper we discuss how to modify the Gadget-2 and evaluate how 
reliable this modified code is, performing a series of tests. 

As we have shown in the previous sections, our tests are successful: numerical errors 
from our changes on Gadget-2 are neglected and we can use this modified code to 
probe alternative theories of gravitation. For example, we can use arguments based 
on galactic dynamics to exclude or not the Yukawa potential hypothesis. This is very 
important, due to the fact that alternative theories appear frequently in the literature, 
but so far there are no many investigations in the galactic scales using JV-body codes 
and galactic dynamics formalism. We will show, in future papers, how to test these 
alternative theories using early and late-type galaxies. 

Acknowledgements The authors would like to thank the Brazilian agencies CAPES, CNPq 
and FAPESP for support. The authors would like also to thank Marcio Eduardo da Silva Alves 
for interesting suggestions, discussions and for the critical reading of the paper. 

References 

1. Piazza, F. & Marinoni, C, Physical Review Letters, 91, 14 (2003). 

2. Signore, R. L., Monthly Notices of the Royal Astronomical Society, 364, 1219 (2005). 



17 



3. dc Araujo, J. C. N., & Miranda O. D., General Relativity and Gravitation, 39, 777 (2007). 

4. Toomre, A., The Astrophysical Journal, 138, 385 (1963). 

5. Rodriguez-Meza, M. A., Cervantes-Cota, J. L., Pedraza, M. I., Tlapanco, J. F., de la 
Calleja, E.M., General Relativity and Gravitation, 37, 823 (2005). 

6. Cervantes-Cota, J. L., Rodn'gucz-Meza, M. A., Gabbasov, R., Klapp, J., Revista Mexicana 
dc Fisica, S53, 22 (2007) 

7. Rodn'gucz-Meza, M. A., Gonzalez-Morales, A. X., Cervantes-Cota, J. L., Gabbasov, R., 
Journal of Physics: Conference Series, 91, 012012 (2007) 

8. Brans, C, Dicke, R., Physical Review, D124, 925 (1961). 

9. Springel, V., Monthly Notices of the Royal Astronomical Society, 364, 1105, 4 (2005). 

10. Visser, M., General Relativity and Gravitation, 30, 12 (1998). 

11. Barnes, J. E. & Hut, P., Nature, 324, 446 (1986). 

12. Goldstein, H.: Classical Mechanics. Addison- Wesley Publishing Company (1980). 

13. Springel, V. & White, S. D. M., Monthly Notices of the Royal Astronomical Society, 307, 
162 (1999). 

14. Binney, J., Tremaine, S.: Galactic Dynamics, Princeton: Princeton Univ. Press (2008) 

15. Gabbasov, R. F., Rodriguez-Meza, M. A., Klapp, J., Cervantes-Cota, J. L., Astronomy 
and Astrophysics, 449, 1043 (2006). 



